Bidirectional transport in a multispecies TASEP model 
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We study a minimal lattice model which describes bidirectional transport of "particles" driven 
along a one dimensional track, as is observed in microtubule based, motor protein driven bidirectional 
transport of cargo vesicles, lipid bodies and organelles such as mitochondria. This minimal model, a 
multi-species totally asymmetric exclusion process (TASEP) with directional switching, can provide 
a framework for understanding the interplay between the switching dynamics of individual particles 
and the collective movement of particles in 1-dimension. When switching is much faster than 
translocation, the steady state density and current profiles of the particles are homogeneous in the 
bulk and are well described by a Mean-Field (MF) theory, as determined by comparison to a Monte 
Carlo simulation. In this limit, we construct a non-equilibrium phase diagram. Away from this 
fast switching regime, the MF theory fails, although the average bulk density profile still remains 
homogeneous. We study the steady state behaviour as a function of the ratio of the translocation 
and net switching rates, Q, and find a unique first-order phase transition at a finite Q associated 
with a discontinuous change of the bulk density. When the switching rate is decreased further 
(keeping translocation rate fixed), the system approaches a jammed phase with a net current that 
tends to zero as J ~ 1/Q- 



I. INTRODUCTION 



Intracellular transport of molecules or cargo vesicles 
often proceed bidirectionally along one-dimensional cy- 
toskeletal filaments such as microtubules. The most well 
studied of these are the dynamics of pigment granules in 
melanophores [l| , mitochondria in the axons of neuronal 
cells , endosomes in Dictyostelium cells [|[ and lipid- 
droplets in Drosophila embryos 0, 0. Bidirectionality is 
achieved by binding to oppositely directed motor pro- 
teins such as kinesin and dynein working in a regulated 
manner 0, H| . Typically intracellular transport along a 
single filament involves many cargo vesicles /organelles, 
having a distribution of sizes. In order to have unim- 
peded and efficient movement of the cargo, the "single- 
particle" switching mechanism must be coordinated with 
the collective "many-particle" behaviour. This provides 
the motivation for the study of minimal models of bidi- 
rectional transport that highlight this interplay between 
individual switching and collective movement of motile 
elements. 

With this in mind, we extend the popular one- 
dimensional (Id) driven diffusive models to describe bidi- 
rectional movement. In spite of their inherent simplic- 
ity, driven diffusive systems exhibit a rich phase diagram 
characterized by novel non-equilibrium steady states with 
finite macroscopic current [9l-ll2| and boundary induced 



phase transitions 0, [IS- 15]. A particular example of 
such models is the Totally Asymmetric Exclusion Pro- 
cess (TASEP), in which the 'particles' hop unidirection- 
ally with a single rate and interact via hard-core re- 
pulsion. Introduction of oppositely directed species in 
the TASEP model would inevitably lead to jamming. 
A natural way to avoid this would be to have stochas- 
tic (or regulated) directional switching of the oppositely 
directed species. Here we study a two species TASEP 
model for driven transport where the two species are op- 
positely directed and can interconvert stochastically lead- 
ing to directional switching. This minimal model shows 
bidirectional transport and a variety of nonequilibrium 
phase transitions. We focus our attention on the situa- 
tion where the translocating particles remain attached to 
the 1-d filament, the conserved model, first introduced 
in [l6j], in contrast to the situation studied in [lj], fl8| . 
where translocating particles are subject to Langmuir 
(un)binding kinetics leading to particle nonconservation 
in the bulk. 

In Section [HI we describe the model and study it us- 
ing a continuum mean-field (MF) approach in the limit 
for which the dynamics of switching is much faster than 
the translocation of the individual species. A compar- 
ison of the mean field results with Monte- Carlo simu- 
lations show very good agreement for the density and 
current profile in the bulk. We also find that the steady 
state solution for average density and current is always 
homogeneous. In Section ILUl we use a domain wall anal- 
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ysis [13| to study perturbations of the steady states and 
construct a nonequilibrium phase diagram. The diverse 
nature of phases are analyzed in terms of variations of 
switching rate. The model also exhibits switching rate 
induced current reversal. Away from this fast switching 
regime, MF theory fails. In Section IPVI we use numerical 
simulations to systematically study departures from this 
fast switching limit. We study it in terms of the ratio of 
the translocation and net switching rates, Q. We find a 
first-order phase transition, associated with a discontinu- 
ous change of bulk density on variation of Q. For a fixed 
translocation rate when Q is increased further, the sys- 
tem approaches a jammed phase with a net current that 
approaches zero as J ~ 1/Q and a density profile that is 
independent of the boundary (un)loading rates. Finally 
in Section|V]we discuss the possible ramifications of these 
results on filament based intracellular transport. 



various dynamic processes and the corresponding rates. 
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II. TWO SPECIES MODEL 



FIG. 1: Dynamical rules for the minimal model involv- 
ing translocation, switching and entry/exit. The reversible 
switching rates between (+) and ( — ) are also shown. 



Cargo vesicles bound to multiple +/— end directed 
motors moving along a microtubule, either make steps 
towards the +/— end or stop, as a result of the regu- 
lated switching of the bound motors. Thus, transitions 
in the cargo vesicle through 'chemical space' (defined by 
the states of the multiple motors) gives rise to movement 
in physical space (a translocation towards the +/— end 
or stationary). For a minimal model, we will restrict the 
instantaneous state of a 'particle' to be either in a (+)- 
directed state or (— )-directed state, with fixed switching 
rates between them [l9[ . Switching can be a result of an 
internal regulatory mechanism attributed to self regula- 
tion by motors or external regulation [3|. The particle 
moves along the 1-dimensional filament of length L, rep- 
resented by a lattice of N sites, labeled i — 0,. . . ,N — 1, 
with lattice spacing e = L/N. The sites i — and 
i = N — 1 define the left and the right boundary, re- 
spectively. Each site i is occupied by either (+), ( — ) or 
a vacancy (0) . 

The (+)/(—) hop to the (right/left) with the same rate 
a, if and only if the adjacent site is vacant. The rate 
of switching /c+(fc_) corresponds to interconversion rate 
from a right (left) moving state to left (moving) moving 
states. The (+) enters the left boundary with a rate cn+ 
if that site is vacant, and leaves the right boundary with 
a rate f3 + . The (— ) species enter the right boundary 
with a rate a_ if that site is vacant, and leaves the left 
boundary with a rate /3_ . Figure 1 is a schematic of the 



The discrete equations motion for the left and right 
moving vesicles at site i are expressed in terms of occu- 
pation number operators, 



dt 



for i 



i(l 



i±l 



-fcfcl)] 
(1) 



1, . . . , N — 2. The terms on the right hand side 
have their usual interpretation of gain and loss terms 
at each site i, arising from translocation and switching 
processes. The corresponding current for the left and 
right moving particles can be expressed as 
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The vacancy occupation number nf is determined by the 



o 



conservation law n\ 



The boundary conditions can be expressed in terms 
of the instantaneous vacancy number, n?. At the left 
boundary site, i = 0, a vacancy enters, when either a (— ) 
particle leaves the left boundary or a (+) particle hops 
to the neighbouring site in the bulk. Similarly a vacancy 
leaves the boundary site at the left boundary when either 
a (+) enters the boundary site, i = or a (— ) particle 
enters this site the next neighbouring site in the bulk, 
i = 1. Accordingly, 



dn° 

—j-j- = ari^ n 1 + P-n Q — a 



+n — an 1 n 



(3) 
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Similarly the evolution of vacancies at the right boundary 
reads 



dn%_ 1 



dt 



= an+_ 1 n%_ 2 +f3 + n+_ 1 -a-n%_ 1 -an+_ 2 n° N _^ 

(4) 



If one averages the occupation number of each species 
over the events that may occur between t and t + dt and 
over all histories up to time t, then the corresponding 
equations of motion for the expectation value of the oc- 
cupation numbers are given by, 

= a[(A Tl ^)-{nt^)]-k ± {nt) + k T { ni \$) 

For the boundary at i = 0, 

d{<) 



dt 



= a i n Q n i) + /M«o > - a +( n o) - a ( n i n o) ( 6 ) 



For the boundary at i = N — 1, 



dt 



o{^-l n W-2) +0+( n N-l) - a -( n N-l) 



( n N-2 n N-l) 



(7) 



To determine the expectation value of occupation number 
one requires the knowledge of all higher order correlations 
which makes it analytically intractable jllj . In specific 
cases like TASEP, owing to special conservation laws, this 
moment hierarchy can be handled exactly using matrix 
methods [2(| [2l|- In situations where the fluctuations 
of the number densities are smaller than mean density, 
a mean field approach works remarkably well [TP . fl3j . 
We will see that mean field approach works very well 
when the process of interconversion is much faster than 
the translational hopping on the lattice, i.e., fc+, fc_ 3> a 
(referred to as the fast switching limit). However away 
from this limit, the mean field approach fails, as seen by 
comparing with results from Monte-Carlo simulations. 



Mean field and Continuum limit 



Defining p 



+ — „ + 



and p: 



the mean field approx- 



imation amounts to factorizing all the two-point correla- 
tors arising out of the different combinations of , n~ 
as a product of their averages, 



(nfnf + 



i) = ( 



nf)( n t+i) = pfpf+i 



(8) 



After taking the expectation values of the operator equa- 
tions, we go over to a continuum description. We normal- 
ize the total length of the lattice, L to 1 and let N — > oo 
so that the lattice spacing e = — > in the thermody- 
namic limit. The rescaled position variables are defined 
as x = i/N — 1, < x < lTogo from discrete number 
representation to continuum density, we Taylor expand, 



(nf)= p±{x) 
nf +1 ) 



P±(x) + ep±(x) + -^p±{x) + . 



p ± (x) - ep±(x) + —p±{x) + . 



(9) 



It useful to re-express the two species densities p+(x) 
and p- (x) in terms of the total cargo density and the 
relative concentration of the different species, p = p++p- 
and 4> = p+ — p-. To leading order in e, the evolution 
equations then read as, 



(10) 



— = -k + (p + ct>) + k-{p-(j))-ea— [p(l-p)fll) 



At steady state, Eq.(fT0| implies that 0(1 — p) is con- 
stant in space. We note that — ap{l — p) is the expression 
for vacancy current. At steady state, in the absence of 
any sink and source terms in the bulk, the conservation 
of vacancies in the bulk implies that the vacancy current 
is spatially homogeneous. This implies that both 4> & n d 
p are individually spatially homogeneous at steady state. 

In the thermodynamic limit of e — > 0, following 
Eq. flTTj) . the two densities become proportional to each 
other. 



(k- - k+)p _1-Kj i 
(k+ + k-) ~ W+ 



■p{x) 



(12) 



where we have introduced the effective ratio of intercon- 
version rates K + = r^- and W+ = 1 + K + . From this 
relation we obtain the corresponding relation between the 
local occupancy of each species 



p + (x) = K + p-(x) 



(13) 



Using eq. (|13p . the corresponding expressions for the 
steady state currents for the individual species are, 



J+ = a p+(l - W+p+) 
J- = -ap-(l-W+p+) 



(14) 
(15) 



4 




- (b) 


1 1 1 1 1 1 
• Simulation (+) 
■ Simulation (-) 
— ■ Mean Field (+) — 
Mean Field (-) 






. 1 


1,1,1," 



0.2 0.4 0.6 0.1 
X 



1 0.2 0.4 0.6 0.1 
X 



0.4 
0.3 
-0.2 
0.1 


-0.1 



1 



This also follows as a consequence of Eq.© and Eq.® af- 
ter appropriately invoking the Mean-Field limit and then 
taking the continuum limit. At steady state R en t = Rext- 
Using the fact that p + = K + p_ we obtain, 



P+(0) 



M =F [M 2 - Aa+W+ (1 -K+)\ 
2W+ (1 - K+) 



(18) 



where M = a+W+ + f3-K+ + 1-K+. Similarly, at the 
right boundary of the filament, x = 1, 

p+{1) = r+[* + 4"f+lX-K+)]* (19) 

H+y 1 2W+(1-K+) v ' 

where P = -a„W+ -0 + + l- K+. 



FIG. 2: Comparison of the steady state density and current 
profiles obtained by Monte Carlo with the Mean Field theory 
in the fast switching limit, shows very good agreement be- 
tween Monte- Carlo simulation and MF results. These steady 
state profiles are for the High Density (HD) phase, (a) Profile 
for Density of (+) and (-) (b) Normalized current: J = J±/a. 
Here a+/a = 0.3, /3+/a = 0.2, a_/a = 0.4, /3-/a = 0.1 and 
K+ =0.5. Monte Carlo simulation is done for a system size 
N = 500. 



Eq. fTU)) also implies that the currents J+, J_ are con- 
stant in the bulk with no spatial variation. Further, they 
are related to each other by, 

J_ = -K+J+ 

The total particle current in the bulk reduces to, 

J(x) = J+(x) + J_(x) = (1 - K+)J+{x) (16) 

which explicitly shows that switching dynamics due to 
regulatory mechanism allows for current reversal. For 
< K + < 1, the overall flux of vesicles is from left to 
right while for K + > 1, the overall flux is from right to 
left. 



III. PHASES AND DOMAIN WALL ANALYSIS 



In the limit of e — > 0, both the boundary conditions 
cannot be satisfied simultaneously. Therefore either the 
steady state solution satisfies the boundary condition at 
the left (x = 0) or it satisfies the boundary condition at 
the right (x = 1) and both of these correspond to the 
two possible stationary states and correspondingly two 
different phases. 

In addition, there is a one more phase, the maximal 
current phase. It corresponds to a situation where the 
stationary state is independent of the entry and exit rates 
of particles at the boundaries and the state is determined 
by maximizing the current in the lattice. Maximizing the 
expression for currents in eq. (114|) , the density and current 
in this phase are, 



1 



f>+ 



J+ = 



2W+ 
1 

4VT+ 



(20) 



B. Boundary Condition 



The steady state condition at the boundaries, is that 
the entry and exit rates of vacancies should match. In 
the thermodynamic limit, the entry and exit rates of the 
vacancies at site i = read, 

R ent = p + (0)[l-p + (0)-p_(0)]+/3_p_(0) 
Re*t = [a++p-(p)][l-p+(p)-p-(p)] (17) 



Thus, similar to TASEP, there are three possible 
phases : Low density (LD) phase corresponding to sta- 
tionary state satisfying the left boundary condition at 
x = 0, a High density (HD) phase corresponding to sta- 
tionary state satisfying the right boundary condition at 
x = 1 and a maximal current (MC) phase correspond- 
ing to the stationary state which is independent of the 
boundary conditions. 
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FIG. 3: Phase diagram in the a+ — P+ plane, for fixed values 
K+ = 0.5, P- = 0.1. (a) q_ = 0.4 and (b) a_ = 0.6. Beyond 
the multicritical point at a*_ = 0.5, there is no MC phase, 
like in case of (b). 



the equation of the MC-HD phase boundary, 



(23) 



which is a straight line parallel to j3+ axis. We see that 
beyond a critical value of a_ (a multicritical point ), the 
maximal current phase does not occur in physical range 
of entry and exit rates. Thus the high density phase 
and the low density phase are the only two phases which 
can occur in such a situation. Finally, the LD-HD phase 
boundary is given by, 



P+ = (1 - K+ - a_W+) + 



4q_(1-^)-Z 2 
2Z 



(24) 



where, Z = 2[^- - p+(0)](l - K%). Note that unlike 
TASEP, the phase boundary separating the low and high 
density phases is not a straight line. Fig. 3 shows the 
mean field phase diagrams and effect of control parameter 
in expelling the maximal current phase. 



In order to construct the phase diagram, we need a 
selection criterion for a given phase. A reasonable se- 
lection criterion is provided by the domain wall analysis 
[l3|. A domain wall is initially constructed as an inter- 
facial region separating two distinct stationary states. If 
this configuration is unstable, then the domain wall will 
move to totally eliminate the unstable state. At phase 
coexistence, the domain wall velocity will be zero. 

The expression for the velocity of the domain wall is 

J5-Ji 



v d 



(21) 



where J+,p+ are the density and current to the right of 
the domain wall and , p\ are the density and current 
to the left of the domain wall. 



To obtain the LD-MC phase boundary, we compute 
the velocity of the LD-MC domain wall. Taking p~t = 

Pmc an d p\ = Pld ' an d putting Vd = in eq. (|2"Tj) . we 
obtain the equation for the LD-MC phase boundary, 



2 Wa 



(2/3_K+ +1-K+ 



(22) 



which in the a + — /3+ plane, is a straight line parallel 
to a + axis. To obtain the MC-HD phase boundary, we 
compute the velocity of the MC-HD domain wall. Taking 
Pl = P~mc an d P~r = Phd^ ano - putting Vd = 0, we obtain 



IV. BREAK DOWN OF MEAN FIELD AND 
APPROACH TO JAMMING 



So far our study has been restricted to the fast switch- 
ing regime, where the rates of interconversion are faster 
than translocation. In order to study deviations away 
from this particular limit, we define a dimensionless pa- 
rameter Q = k+ ^_ k , which is the ratio of the transloca- 
tion rate and the net switching rate. Thus fast switching 
corresponds to the limit of Q = 0. 

Most realistic cellular situations are far from this fast- 
switching limit. To see this, we estimate the typical 
values of Q for transport of lipid-droplets in wild type 
Drosophila embryos [6| . The average run velocity of the 
lipid droplet, before a switching or pause event, is mea- 
sured to be ~ 0.2/ims -1 for short runs (0.03 — 0.1 /im) 
and ~ 0.5/ims -1 for long runs (0.5 — 1 fim) Q. From 
this, we estimate the switching rates to be in the range 
of ~ (1 — 4)s _1 . Assuming a step size for a single kinesin 
motor be 8nm and using the measured values for average 
run velocity, we estimate the typical translocation rate a 
to be ~ (25 — 50)s _1 , This would imply that the typical 
range of Q would be between 3 — 25. 

Can we employ the same theoretical mean field strat- 
egy, as in the Q = limit ? Unfortunately we find that 
the mean field predictions fail completely when Q > 0, as 
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FIG. 4: (a) Density profile of p + as a function of Q — fc+ ° fc , 
the ratio of translocation and switching rate, which measures 
the deviation from the fast switching limit. The dashed line 
is the MF prediction for density in HD phase, (b) Plot of the 
corresponding discontinuous change of bulk density at Q = 
0.84, a first order transition from a LD phase to a HD phase. 
Here a+/a = 0.2, j3+/a = 0.5, a_/a = 0.4, /3_/a = 0.8 and 
K+ — 0.25. Monte Carlo simulation is done for a system size 
N = 1000. 



FIG. 5: (a) Plot of vacancy density as a function of Q is 
shown, (b) Variation of the current of the (+) species, J+ 
as a function of Q is shown. The current decreases with in- 
creasing Q, diverging strongly from mean field predictions 
and asymptotically approaching zero for Q — > oo . Inset fig- 
ure shows the plot of log Q vs log J+ for the same data as 
in (b). A straight line fit to the inset plot for data points 
for which Q > 2 shows that J ~ Q", where v ~ —1.097. 
Here a+/a = 0.3, f}+/a = 0.5, a_/a = 0.4, /3_/a = 0.1 and 
K + — 0.5. Monte Carlo simulation is done for a system size 
N = 400. 



is revealed by comparison with Monte Carlo simulations. 
Our 'exact' analytic study in the Q = limit is never- 
theless useful in providing a conceptual basis to analyze 
the numerical simulations when Q > 0. 

In Fig. 4 we show the density profile obtained for dif- 
ferent values of Q, keeping the entry and exit rates of 
particles at the boundaries fixed. Fig. 4(a) shows how 
for values of Q close 0, the MF prediction for densities 
matches well with the Monte-Carlo simulations. The pa- 
rameter values shown in the figure correspond to a region 
in Low Density (LD) phase. On increasing Q, the bulk 
density increases continuously and starts deviating from 
the MF prediction. However the ratio of densities of the 
two species expressed in eq.(fT5]) remains unaltered with 
change in Q. The corresponding value of current in the 
system also deviates from the MF value and decreases 
continuously with the increase in Q. Fig. 4(b) shows that 
at Q rj 0.84, there is a first order transition, associated 
with a discontinuous change of density from a low den- 
sity (LD) phase to a high density (HD) phase. The cur- 
rent remains continuous across the phase transition. This 
non-equilibrium phase transition is unique, arising out of 
an interplay between translocation and switching dynam- 
ics in this driven system. However if the initial phase at 



Q = is cither a HD phase or maximal current (MC) 
phase then no such discontinuous change of density is 
observed on increasing Q. Instead the bulk density in- 
creases continuously with increase in Q, saturating to a 
maximum value. 

Fig. 4(a) shows how the value of the density in the HD 
phase is significantly higher than the MF prediction. On 
further increasing Q, the density of (+) species increases 
and asymptotically approaches a maximum value. Cor- 
respondingly as shown in Fig. 5(a), the vacancy density 
reduces. With an increase of Q, the net current in the 
system approaches zero. Further, the steady state also 
becomes independent of the boundaries for higher values 
of Q. From the log- log plot it can inferred that the net 
current goes as, J ~ l/Q, for large Q. This would im- 
ply that J ~ k±; so that when the switching rates are 
much smaller than the translocation rates, then the net 
particle flux in the lattice is controlled by the switching 
process alone. 

However, as far as we can tell from our numerical stud- 
ies of this model, we do not observe any jamming phase 
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transition at finite Q, upto the values of Q = 500 that we 
have investigated. Beyond that due to the smallness of 
values of the current, statistically reliable values of cur- 
rent is difficult to obtain numerically. Figure 5(b) shows 
how J + reduces and system approaches jamming. Al- 
though the average steady state profile is homogeneous, 
there is an increased tendency for the vacancies to clus- 
ter, as the system approaches jamming. 

It is interesting to compare the steady state properties 
of the current model with other closely related model. In 
particular, a two species model proposed in [22j . with dy- 
namic rules : +0 — !• 0+, > —0 , +— ^ — h, has some 

similarity with the dynamics of the model discussed here. 
It was found that for a closed ring, for a certain range of 
relative rates for the process +— ^ — +, there is a phase 
transition associated with the formation of three distinct 
domain which either move with constant drift velocity or 
remain in a pinned configuration, depending on the ra- 
tio of two species of particles. Apart from the nature of 
the boundary conditions, the crucial difference between 
the two models is the aspect of Charge-Parity (CP) con- 
servation - while charge-parity is conserved in |22j |. it 
is explicitly broken in our case because of the switching 
dynamics. We would like to explore whether the model 
that we discuss in this paper can also support such novel 
steady states, if we change the boundary conditions from 
open to closed. 

Based only on our current numerical analysis it is dif- 
ficult to assess whether this model has a finite Q phase 
transition to a jammed state with zero particle current or 
simply a crossover to jamming - clearly this issue needs 
to addressed analytically. 



V. CONCLUSIONS 



In this paper, we study a minimal lattice model which 
describes bidirectional transport of "particles" driven 
along a one dimensional track, as is observed in micro- 
tubule based, motor protein driven bidirectional trans- 
port of cargo vesicles, lipid bodies and organelles such as 
mitochondria. This model can be described as a multi- 



species totally asymmetric exclusion process (TASEP) 
with directional switching and which obeys overall parti- 
cle number conservation in the bulk. It provides a simple 
framework for understanding the interplay between the 
switching dynamics of individual particles and the col- 
lective movement of particles in 1-dimension. 

The model exhibits various interesting features such 
as current reversal of the net cargo, depending on the 
effective ratio of the switching rate K + . This aspect of 
the model is consistent with flux reversals observed in the 
transport of pigment granules in melanophore cellular ex- 
tracts, wherein changing the hormonal conditions leads 
to alterations in the (un)binding rates of the motors 
The model predicts a multicritical behaviour; the maxi- 
mal current phase is absent for the entire physical range 
of (a + — j3 + ) phase plane beyond a critical value of a_. 
Switching rate induced phase transition is another inter- 
esting feature exhibited by this model : thus changes in 
the translocation rate or switching rates of the cargo us- 
ing say biochemical methods, can lead to sharp changes 
in the average density profile of transported cargo. 

A possible extension of the present model is to go be- 
yond local particle conservation by explicitly allowing 
particles to exchange with the environment. Such a sit- 
uation would typically arise when cargo vesicles moving 
along microtubules, detach from it and move along the 
surrounding actin meshwork for a short run time, before 
reattaching to the same microtubule at a different lo- 
cation. Since the individual actin filaments in the actin 
meshwork are short and orientationally uncorrelated, the 
movement of the cargo vesicles in the actin meshwork is 
effectively diffusive. This suggests that the dynamics of 
cargo vesicles along a microtubule embedded in a con- 
fined isotropic actin meshwork can be modeled by an 
appropriate non-conserved version of the current model. 
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